# load("xxx.RData")



### PRODUCTION FUNCTION

library(party)

attach(articles)

# Regression tree with average prior impact:

temp <- ctree(impact ~ 
	avg_impact +
	apmic + methy + macro + bzfin + agloc + num_areas +
	tot_proj + team_deg +
	num_authors + skill_deficit,
	controls=ctree_control(mincriterion=0.99) )

# Plot (Fig. A4):

png("reg_tree_full.png", width=2000, height=500)
plot(temp, type="simple",
  inner_panel=node_inner(temp, digits=2, fill=rgb(0, 0, 0, alpha = 0.02), id=F, pval=F),
  edge_panel=edge_simple(temp, digits=2),
  drop_terminal=FALSE,
  terminal_panel=node_terminal(temp,digits=2, fill=rgb(1, 1, 1, alpha = 0.75), id=T)
  )
dev.off()

# Residualize from effect of average prior impact:

temp <- ctree(impact ~ 
	avg_impact,
	controls=ctree_control(mincriterion=0.95) )

impact_res <- impact - as.numeric( Predict(temp) )

# Regression tree estimate of production function:

prodfn <- ctree(impact_res ~ 
	apmic + methy + macro + bzfin + agloc + num_areas +
	tot_proj + team_deg +
	num_authors + skill_deficit,
	controls=ctree_control(mincriterion=0.99) )

# Plot (Fig. 2):

png("reg_tree.png", width=750, height=500)
plot(prodfn, type="simple",
  inner_panel=node_inner(temp, digits=2, fill=rgb(0, 0, 0, alpha = 0.02), id=F, pval=F),
  edge_panel=edge_simple(temp, digits=2),
  drop_terminal=FALSE,
  terminal_panel=node_terminal(temp,digits=2, fill=rgb(1, 1, 1, alpha = 0.75), id=T)
  )
dev.off()

# Clean up:

detach(articles)

rm(temp)
